Read data
source("scripts/myRLib.R")
funcset <- readLines("output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/func_subset.set")
funcset <- parseFuncSet(funcset)
funcset.names <- names(funcset)
xtyi.files <- strsplit("output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p1.0000e+00.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p3.0000e-01.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p1.0000e-01.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p3.0000e-02.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p1.0000e-02.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p3.0000e-03.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p1.0000e-03.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p3.0000e-04.assoc.logistic:output-direction_XT_YI/fastInsulin__wtccc_t2d__wtccc_ctrl/xtyi_LDpred_p1.0000e-04.assoc.logistic",
":")[[1]]
ps <- strsplit("1,0.3,0.1,0.03,0.01,0.003,0.001,0.0003,0.0001",
",")[[1]]
inter.files <- strsplit("output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p1.0000e+00.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p3.0000e-01.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p1.0000e-01.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p3.0000e-02.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p1.0000e-02.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p3.0000e-03.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p1.0000e-03.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p3.0000e-04.assoc.logistic:output-run_step2/fastInsulin__wtccc_t2d__wtccc_ctrl/interaction.func_subset_LDpred_p1.0000e-04.assoc.logistic",
":")[[1]]
prs.files <- strsplit("../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p1.0000e+00.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p3.0000e-01.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p1.0000e-01.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p3.0000e-02.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p1.0000e-02.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p3.0000e-03.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p1.0000e-03.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p3.0000e-04.txt:../step1/output-run_step1/fastInsulin__wtccc_t2d__wtccc_ctrl/ldpred_LDpred_p1.0000e-04.txt",
":")[[1]]
Epistasis
for (i in 1:length(xtyi.files)) {
cat("##", funcset.names[i])
cat("
")
cat("
")
snp.subset <- funcset[[funcset.names[i]]]
# read I ~ SNP
prs <- read.table(prs.files[i], header = T)
prs$sid <- as.character(prs$sid)
prs.sub <- prs[prs$sid %in% snp.subset, ]
# read others
xtyi <- read.table(xtyi.files[i], header = T)
xtyi$SNP <- as.character(xtyi$SNP)
yt.xtyi <- xtyi[(xtyi$SNP %in% snp.subset) & (xtyi$SNP %in%
prs.sub$sid), ]
yt.xtyi.xt <- yt.xtyi[yt.xtyi$TEST == "ADD", ]
yt.xtyi.yi <- yt.xtyi[yt.xtyi$TEST == "PRS", ]
inter <- read.table(inter.files[i], header = T)
inter$SNP <- as.character(inter$SNP)
inter <- inter[inter$TEST == "ADDxPRS", ]
yt.inter <- inter[(inter$SNP %in% snp.subset) & (inter$SNP %in%
prs.sub$sid), ]
# order all rows
yt.xtyi.xt <- yt.xtyi.xt[order(yt.xtyi.xt$SNP), ]
yt.xtyi.yi <- yt.xtyi.yi[order(yt.xtyi.yi$SNP), ]
yt.inter <- yt.inter[order(yt.inter$SNP), ]
prs.sub <- prs.sub[order(prs.sub$sid), ]
# check direction
prs.sub$nt1 <- as.character(prs.sub$nt1)
yt.xtyi.xt$A1 <- as.character(yt.xtyi.xt$A1)
yt.inter$A1 <- as.character(yt.inter$A1)
if (sum(prs.sub$nt1 != yt.xtyi.xt$A1) != 0 | sum(prs.sub$nt1 !=
yt.inter$A1) != 0) {
print("Direction does not match. Exit!")
quit()
}
yt.xtyi.xt.mean.std <- getMeanStd(yt.xtyi.xt$SNP, yt.xtyi.xt$P,
yt.xtyi.xt$OR, "Y~Gi+I (Gi)")
yt.xtyi.yi.mean.std <- getMeanStd(yt.xtyi.yi$SNP, yt.xtyi.yi$P,
yt.xtyi.yi$OR, "Y~Gi+I (I)")
yt.inter.mean.std <- getMeanStd(yt.inter$SNP, yt.inter$P,
yt.inter$OR, "Y~Gi*I")
prs.sub.mean.std <- data.frame(SNP = prs.sub$sid, Mean = prs.sub$ldpred_beta,
Std = 0, Type = "I~Gi")
# create data frame for plot
df <- rbind(yt.xtyi.xt.mean.std, yt.xtyi.yi.mean.std, prs.sub.mean.std,
yt.inter.mean.std)
p <- ggplot(df) + geom_point(aes(x = SNP, y = Mean)) + geom_errorbar(aes(x = SNP,
ymin = Mean - 1.96 * Std, ymax = Mean + 1.96 * Std),
colour = "black", width = 0.1) + geom_hline(yintercept = 0,
color = "red") + facet_grid(. ~ Type, scales = "free_x") +
coord_flip() + theme(axis.text.x = element_text(angle = 90,
hjust = 1))
p2 <- ggplot(yt.inter) + geom_bar(aes(x = SNP, y = -log10(P)),
stat = "identity") + coord_flip() + ggtitle("Interaction P-value")
print(p)
print(p2)
cat("
")
cat("
")
}
p1.0000e+00
p3.0000e-01
p1.0000e-01
p3.0000e-02
p1.0000e-02
p3.0000e-03
p1.0000e-03
p3.0000e-04
p1.0000e-04
